% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Function for performing CG bootstrap in BR2020 Application
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function out = bootstrapCG(fwd,pred,initCondFwd,initCondPred,T_all,idxSampleStartPCs,idxSampleEndPCs,p)

        numPred = size(pred,2);
    
        simdretapprox = fwd(1:end-1,2:end) - fwd(2:end,1:end-1);
        varJoint = var1BC([fwd(:,end), pred], p.CG_B_BC,false);        
        Zfwd = varJoint.VBC(2:end,1);
        Zpred = varJoint.VBC(2:end,2:end);
        Zinit = [Zfwd, simdretapprox, Zpred];
        Z = [Zinit; Zinit];
        idxRand = randi(T_all-p.M,[ceil(T_all/p.M),1]);
        idxBootstrap = kron(idxRand,ones(p.M,1)) + kron(ones(ceil(T_all/p.M),1),(0:p.M-1)');
        idxBootstrap = idxBootstrap(1:T_all-1);
        bsZ = Z(idxBootstrap,:);

        bsV = bsZ(:,[1,end-numPred+1:end]);    
        bsV = bsV - ones(T_all-1,1)*mean(bsV);

        tmpJoint = NaN(T_all,numPred+1);       
        tmpJoint(1,:) = [initCondFwd(end), initCondPred];
        for t = 2:T_all
            tmpJoint(t,:) = [1 tmpJoint(t-1,:)]*varJoint.PsiBC + bsV(t-1,:);
        end
    
        out.bsPred = tmpJoint(:,2:end);        
        bsDret = [NaN(1,size(fwd,2)); NaN(T_all-1,1) bsZ(:,2:size(simdretapprox,2)+1)];
        bsFwd = buildForwards(bsDret, initCondFwd, tmpJoint(:,1));
        
        % Make yields, forwards and returns for the Long Sample
        out.bsFwdCG = bsFwd;
        out.bsYldCG = NaN(size(bsFwd));
        out.bsYldCG(:,1) = bsFwd(:,1);
        for i = 2:size(bsFwd,2)
            out.bsYldCG(:,i) = ((p.yldMats(i)-p.yldMats(i-1))*out.bsFwdCG(:,i) + p.yldMats(i-1)*out.bsYldCG(:,i-1))/p.yldMats(i);
        end

        % Make PCs for "Sample"; subtract 4 because we start data in 1960Q4
        [tmp_evecs,tmp_eigs] = eig(cov(out.bsYldCG(idxSampleStartPCs:idxSampleEndPCs,:)));
        [~,idx] = sort(diag(tmp_eigs), 'descend');
        tmp_evecs = tmp_evecs(:,idx);
        bsPCs = out.bsYldCG*tmp_evecs(:,1:3)*diag(sign(tmp_evecs(end,1:3)));
        bsPC1 = bsPCs(:,1);
        bsPC2 = bsPCs(:,2);
        bsPC3 = bsPCs(:,3);
        % Make scaled PCs
        bsSc = [sum(tmp_evecs(:,1)), tmp_evecs(1,2)-tmp_evecs(end,2), tmp_evecs(end,3)-2*tmp_evecs(4,3)+tmp_evecs(1,3)];
        out.bsPCssc = bsPCs*diag(1./bsSc);
        out.bsPC1sc = out.bsPCssc(:,1);
        out.bsPC2sc = out.bsPCssc(:,2);
        out.bsPC3sc = out.bsPCssc(:,3);
                
        [out.bsXrn,out.bsXrn4] = makeReturnsBR(out.bsYldCG);

end